Tensor-entanglement renormalization group approach to 2D quantum systems 
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Traditional mean-field theory is a simple generic approach for understanding various phases. 
But that approach only applies to symmetry breaking states with short-range entanglement. In this 
paper, we describe a generic approach for studying 2D quantum phases with long-range entanglement 
(such as topological phases). Our approach is a variational method that uses tensor product states 
(also known as projected entangled pair states) as trial wave functions. We use a 2D real space RG 
algorithm to evaluate expectation values in these wave functions. We demonstrate our algorithm 
by studying several simple 2D quantum spin models. 



Introduction: To obtain various possible quantum 
phases of a quantum spin system H = ^ -v Si ■ Jij ■ Sj , 
we may use a mean-field approach. The mean-field ap- 
proach can be viewed as a variational approach. For ex- 
ample, to study the possible spin ordered phase of the 
above quantum spin system, we may start with a trial 
wave function Atrial) = ®(«« |T)i + v i where |T); 

and \l)i are spin states on site i. The spin ordered phases 
can be obtained through minimizing average energy by 
changing Ui and t^. But such kind of mean- field the- 
ory only apply to states with short-range entanglement 
(since Atrial) is a direct product state). As a result, we 
cannot use the traditional mean-field theory to under- 
stand quantum phases that have pattern of long-range 
entanglement (such as topologically ordered states and 
other quantum states beyond Landau's symmetry break- 
ing description). [1-3] 

One approach for addressing these phases is to use 
a more general class of trial wave functions known as 
"tensor product states" (TPS) or "projected entangled 
pair states" (PEPS). [4, 5] Tensor product states were 
first discovered in the context of the (ID) density matrix 
renormalization group (DMRG) method [6, 7], but were 
later generalized to higher dimensions and arbitrary lat- 
tices. On the square lattice (Fig. I), the TPS are defined 
by 
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where T^j i is a complex tensor with one physical index 
nii and four inner indices i, j,k,l, ■ ■ ■ . The physical in- 
dex runs over the number of physical states d on each 
site and inner indices runs over D values. Unlike simple 
mean field states, these variational wave functions can 
describe 2D many-body quantum systems [8] with short- 
range entanglement (such as symmetry breaking states) 
as well as long-range entanglement (such as string-net 
condensed states [9]). 

One of main challenges of using this approach in higher 
dimensions is that it is not easy to compute expecta- 
tion values in these states. In this Letter, we describe a 
simple solution to this problem in two dimensions. Our 
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FIG. 1: Tensor-network - a graphic representation of the 
tensor-product wave function (1), (a) on a ID chain or (b) 
on a 2D square lattice. The indices on the links are summed 
over. 



approach - which we call the tensor entanglement renor- 
malization group (TERG) method - is an approximation 
scheme based on the 2D real space RG method devel- 
oped in Ref. 10. A different real space RG method can 
be found in Ref. [11]. 

As we mentioned earlier, this kind of variational ap- 
proach has the advantage that it can potentially address 
2D quantum many-body states that contain both symme- 
try breaking orders and topological orders. In this paper, 
we will just introduce our algorithm by studying a few 
simple quantum models and compare our results with 
those obtained through other previous methods. The 
application of TERG approach to topologically ordered 
states will be presented in future publications. 

To see the efficiency of the TERG method, let us com- 
pare it with other variational methods for 2D gapped 
systems (see Appendix for an explanation): 



Method 


Error 


VQMC 

ID approach 

TERG 


e-l/T 1 ^ 

e ~ exp(— const ■ logT) 
e ~ exp(— const • (logT) 2 ) 



Here T the calculation time and e is the achieved accu- 
racy of the calculated average energy for a given many- 
body state. The acronym VQMC stands for variational 
Quantum Monte Carlo, while "ID approach" refers to an 
approximation scheme where one replaces the infinite 2D 
lattice by an L x oo lattice, for L large but finite, and 
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FIG. 2: (a) The graphic representation of the inner product 
of TPS, (*|*), in term of the tensor T or the double-tensor T 
(see (2)). (b) After summing over m and identifying (i,i r ) —* 
a i Cm') ~ * Pi we obtain the double-tensor T from T and T* . 



then computes expectation values using a transfer matrix 
approach. 

The TERG method: Let us consider a system with 
translationally invariant H = J^i^i- Hi can always be 
expressed as a summation of products of local operators: 
Hi = 0° + 0\&] + • • • • A key step in the variational 
approach is to calculate the norm and the expectation 
value of Hi for the trial wave function $({m,}): 
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= tTr[T<8> T(g>T<g> •••] , 
(V\Hi\*) = tTr [T° (g> T (g) T (g> • • • ] 

+ tTr [T] <g> Tf <g> T ® •••]+•• ■ 



(2) 



where double-tensors T, T°, a = 0,1,2, are defined as 
(see Fig. 2b) 
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with O'L , the matrix elements of local operators O a in 
the local basis |m). The tensor-trace (tTr) here means 
summing over all indices on the connected links of tensor- 
network (see Fig. 3a). Note that the inner product is 
obtained from a uniform tensor-network. The average of 
the on-site interaction is obtained from a tensor-network 
with one "impurity" tensor T° at site i (while other site 
has T). Similarly, for two-body interactions, the tensor- 
network has two "impurity" tensors T 1 and T 2 at i and 
3- 

Calculating the tensor-trace tTr is an exponentially 
hard calculation in 2D or higher dimensions. Motivated 
by the tensor-renormalization approach developed in Ref. 
[10], we can accelerate the calculation exponentially if 
we are willing to make an approximation. The basic 
idea is quite simple and is illustrated in Fig. 3. Af- 
ter finding the reduced double-tensor T", we can express 
tTr[T<g>T • • • ] ps tTr[T"®T" • • • ] where the second tensor- 
trace only contain a quarter of the double-tensors in the 
first tensor-trace. We may repeat the procedure until 
there are only a few double-tensor in the tensor-trace. 
This allows us to reduce the exponential long calculation 
to a polynomial long calculation. 
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FIG. 3: The indices of the double-tensor have a range D . 
After combine the two legs on each side into a single leg, the 
four linked double-tensors in (a) can be viewed as a single 
double-tensor T' whose indices have a range D 4 . (c) T' can 
be approximately reduced to a "smaller" double-tensor T" 
whose indices have a range D 2 and satisfies tTr[T' ® T' • • • ] ps 
tTr[T" ® T" ■ ■ ■]. 




FIG. 4: (Color online) (a) We represent the original rank-four 
tensor by two rank-three tensors, which is an approximate de- 
composition, (b) Summing over the indices around the square 
produces a single tensor T'. This step is exact. 



The detail implementation of the above TERG ap- 
proach is actually a little more involved. For an uniform 
tensor-network Fig 3a, we can coarse grain it in two steps. 
The first step is decomposing the rank-four tensor into 
two rank-three tensors. We do it in two different ways on 
the sublattice purple and green (see Fig. 4a). On pur- 
ple sublattice, we have T af3lll/ = J2y S\nwy'Szap-y> and 
on green sublattice, we have = J2<y> Siva-i 1 Sif3n-y' ■ 

Note that a, ft, /i, v run over D 2 values while 7' run over 
D 4 values. 

Next we try to reduce the range of 7' through an 
approximation. [10] Say, on purple sublattice, we view 
as a matrix M^.^ v = T Qj3(J „ and do singular 
value decomposition M rcd = UKV^ . We then keep only 
the largest D cut singular values A 7 and define Si^ vl = 

V^r^y uvi ^3ap-y = \Z^yU a f3 tl . Thus, we can approxi- 
mately express T a/ 3 Mt/ by two rank-three tensors S±, S3 



T Q /3 M y — ^ ] *S'3a/37*S'l/jiy7- 
7=1 



(4) 



Similarly, on green sublattice we may also define T a ^„ as 
a matrix and do singular value decompositions, 

keep the largest D cut singular values and approximately 
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FIG. 5: (Color online) Iteration of tensor-network produces a 
coarse grained tensor-network. 



express T a ^ M „ by two rank three tensors S?, S4. 



D cut 
7=1 



(5) 



After such decompositions, the square lattice is de- 
formed into the form in Fig. 4b (see also Fig. 5). The 
second step is simply contract the square and get a new 
tensor on the coarse grained lattice. 

^T-erAp — J]] Sip al S2 ll j3 rJ S3, ytl xSi al ,p (6) 
aft [iv 

The range of indices for the reduced double-tensor T' is 
only £> cu t which can be chosen to be D 2 or some other 
values. Repeat the above two steps twice, we can get the 
reduction from Fig. 3a to 3c. 

The above TERG procedure can be easily generalized 
to tensor-network with "impurities", such as the one in 
Fig. 5 which has four "impurity" tensors. Evaluating 
Fig. 5 will allow us to calculate the averages of up to four- 
body nearest-neighbor interactions (which include on-site 
interaction, nearest-neighbor and next-nearest-neighbor 
two-body interactions). TERG procedure is illustrated 
in Fig. 5. We note that the number and the relative 
positions of "impurity" tensors do not change after each 
iteration. So we repeat the same iterative calculation un- 
til there are only a few tensors in the tensor-trace. Thus 
the calculations of the averages of local operators is also 
reduced to polynomial long calculations. The total com- 
putational complexity is cost time ~ D® ut \ogN on square 
lattice (N is the total number of sites). For gapped sys- 
tems in the thermodynamic limit, the truncation error 



can be estimated as 6 ~ exp 



-const • (log D cut ) [10] 
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FIG. 6: (Color online) (a) Magnetization along the x direction 
(a x ) versus transverse field h. The derivative of magnetization 
has a singularity around h ~ 3.1, indicating the second order 
phase transition, (b) Magnetization along the z direction (a z ) 
versus transverse field h. In the inset is the log plot of (a z ) 
versus \h — h c \, where h c is the critical field. 



Hi in (2), we can obtain the approximated ground state 
with minimized average energy by adjusting the elements 
in the tensor T. 

Examples: To test our TERG algorithm, we hrst cal- 
culate ground state and its magnetization along x and z 
directions for the transverse field Ising model: 



(7) 



We choose the tensor T in Eq. (1) to be real and has 
90 degree rotational symmetry. We also choose the inner 
dimension D — 2 and keep 18 singular values at each iter- 
ation (D cu t = 18). The total system size is up to 2 9 x 2 9 
sites. The average energy for a tensor T is calculated 
using the TERG approach. We use Powell minimization 
method to find the minimal average energy and the cor- 
responding tensor which gives us the variational ground 
state. 

In Fig. 6 we plot the polarization along x direction 
and z direction in the variational ground state. We note 
that despite the a z — » — a z symmetry in the Hamilto- 
nian, the tensor T that minimize the average energy may 
break the a z — > — <j z symmetry and give rise to non-zero 
polarization in z direction. We find a second order phase 
transition at h c « 3.08. We further fit the critical expo- 
nent 



(a z ) = A\h - h c 



(8) 



After calculating the inner product and the average of 



with j3 w 0.333 ± 0.003. Both the values of critical field 
and critical exponent f} here are very close to the QMC 
results, with h^ MC ~ 3.044[12] and f3 QMC ~ 0.327[13]. 
They are much better than the meanfield results h c = 4 
and = 0.5. 

To see the truncation error caused by D cutl we plot 
the ground state energy (per site) of Eq. 7 as a function 
of D cut , for h = 2.8,3.2, and h = h c = 3.08 (see Fig. 
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FIG. 7: Ground state energies of the transverse Ising model 
for different D C ut- 

7). The energies for different ft's are shifted by different 
constants so that the three curves can be fitted into one 
window. Notice that for off critical systems (h = 2.8, 3.2), 
the energy converges very quickly for small D cut (~ 26). 
Even at the critical point h — h c = 3.08, the error in 
energy per site is of order 10~ 4 . The truncation error is 
much smaller for gapped off-critical states. 

As another more stringent test, we also apply the 
TERG method to study Heisenberg model H = Y^Uj) 
Sj on square lattice which contains gapless excitations. 
Again we choose D = 2, D cut — 18 and total system size 
2 9 x 2 9 sites. We choose tensors Ta and Tb to be real and 
has 90 degree rotational symmetry on sublattice A and 
B. We find the ground state energy to be —0.33 per bond, 
which is quite close to the best QMC results(-0.3350).[14] 
The TERG method also allows us to calculate correlation 
function using tensor-network with two "impurity" ten- 
sors with arbitrary separations. Through the long-range 
correlation function, we find that the total magnetiza- 
tion is m = ^{SfS? + S?SV + S?S?) = 0.39, which is 

larger than the QMC rcsults(0.307).[14] We see that a 
small error in ground state energy (which depends only 
on short-range correlation) can leads to a larger error on 
correlations at long distances. 

Conclusions and discussions: The TERG approach 
is a simple generic method to obtain various quantum 
phases and quantum phase transitions for quantum sys- 
tems in any dimension. The most important feature of 
TERG approach is that it can handle quantum states 
with long-range entanglement (such as topologically or- 
dered states) . When we use traditional mean-field theory 
to calculate quantum phase diagram, the topological or- 
dered phases cannot appear in such a mean-field phase 
diagram, since the mean-field states are limited to those 
with short-range entanglement. The TERG approach 
solves this problem and can generate phase diagrams that 
contain both symmetry breaking states and topologically 
ordered states. 
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Appendix: In a VQMC calculation, the error e is a 
statistical error that scales like 1/N 1 / 2 where N is the 
number of samples. The computational time T scales 
like N. Thus the scaling of the error with computational 
time is given by e ~ l/T 1 / 2 . 

In the ID approach, the error e is a finite size error that 
comes from the truncation of the infinite 2D lattice to an 
L x oo lattice. In a gapped system we expect this error to 
fall off as e~ L l^ where £ is the correlational length. On 
the other hand, the computational time T is exponential 
in L since the method requires diagonalizing a transfer 
matrix whose size is exponentially large in L. We con- 
clude that the error scales with computational time as 

£ ^ — const. -log T 

In the TERG approach, the truncation error for each 
iteration step scales as t\ <~ e -const.-(iogD cut ) ; s i nce cal- 
culating the norm and averages is like calculating the par- 
tition function in Ref. 10. The total truncation error for 
a system of size L is e t ~ (log_L)e~ const -'( logDcut ' since 
such a system requires logL iterations. On the other 
hand, the finite size error is e s ~ e~ L// ^. Minimizing the 
sum of the two errors, we see that the optimal L is given 
by L ~ (log _D cut ) 2 . Since the computational time scales 
polynomially in D cut , we conclude that the total error 
scales like e ~ e -const.-(io g T) ( ne glecting subleading log 
corrections) . 
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